#include <iostream>
#include <string>
#include <fstream>

#include "VortSystem.h"
#include "../tinyxml/tinyxml.h"

using namespace VortPan;
using namespace std;


void WritePoints(std::vector<Vector3D>& Stream, std::string FName){
	ofstream of(FName.c_str());
	
	for (std::vector<Vector3D>::iterator it=Stream.begin(); it!=Stream.end(); it++) {
		of<<it->v[0]<<" "<<it->v[1]<<" "<<it->v[2]<<endl;
	}
	
}

Vector3D ReadVector(TiXmlElement *pEl){
	Vector3D v=V0();
	
	v.v[0]=atof(pEl->FirstChildElement("x")->GetText());
	v.v[1]=atof(pEl->FirstChildElement("y")->GetText());
	v.v[2]=atof(pEl->FirstChildElement("z")->GetText());

	if(pEl->FirstChildElement("Normalize")){
		v=Normalize(v);
	}
	
	
	return v;
}



int main (int argc, char * const argv[]) {
	string Name;
	
	if(argc<2){
		Name="case.xml";
	} else {
		Name=argv[1];
	}

	TiXmlDocument Doc;
	assert(Doc.LoadFile(Name.c_str()));

	TiXmlElement *pRoot;
	assert(pRoot=Doc.FirstChildElement());
	
	VortMesh Mesh;
	VortMesh OrigMesh;
	
	double WeldingTolerance=1e-5;
	
	TiXmlElement *pElem;
	if(pElem=pRoot->FirstChildElement("WeldingTolerance")){
		WeldingTolerance=atof(pElem->GetText());
	}
	
	re Scale=1e-0;
	
	TiXmlElement *pMeshes;
	
	assert(pMeshes=pRoot->FirstChildElement("Meshes"));
	
	TiXmlElement *pEl=NULL;
	string MeshFile;
	bool IsWake;
    bool OnlyPos;
    int IndOff;
	
	
	for(pEl=pMeshes->FirstChildElement("Mesh"); pEl!=NULL; pEl=pEl->NextSiblingElement()) {
		VortMesh LoadMesh;

		MeshFile=pEl->GetText();
		Scale=1.0;

		TiXmlElement *pParamElement;

		pParamElement=pEl->FirstChildElement("Scale");
		if (pParamElement!=NULL) {
			Scale=atof(pParamElement->GetText());
		}

		IsWake=false;
		if(pEl->FirstChildElement("Wake")){
			IsWake=true;
		}
		
		LoadMesh.LoadPLY(MeshFile, IsWake, Scale);
				
		OrigMesh.AddMesh(LoadMesh);
	}
	
	for(pEl=pMeshes->FirstChildElement("RawMesh"); pEl!=NULL; pEl=pEl->NextSiblingElement()) {
		VortMesh LoadMesh;
		
		MeshFile=pEl->GetText();
		Scale=1.0;
		
		TiXmlElement *pParamElement;
		
		pParamElement=pEl->FirstChildElement("Scale");
		if (pParamElement!=NULL) {
			Scale=atof(pParamElement->GetText());
		}
		
		IsWake=false;
		if(pEl->FirstChildElement("Wake")){
			IsWake=true;
		}
        
        OnlyPos=false;
        if(pEl->FirstChildElement("OnlyPos")){
			OnlyPos=true;
		}
		

        IndOff=0;
        if(pEl->FirstChildElement("IndOffset")){
			IndOff=atoi(pParamElement->GetText());
		}
		
		LoadMesh.ImportRawData(MeshFile, IsWake, Scale, IndOff, OnlyPos);
		
		OrigMesh.AddMesh(LoadMesh);
	}
	
	
	
	
	OrigMesh.Weld(WeldingTolerance);

	

	

	
	
	
	int NIter=1;
	re RelaxationFactor=1;
	re VortexCoreSize=0.1;
	Vector3D WakeDir=e1();
	re StartSegmentLength;
	re EndDistance;
	re NextElementRatio;

	if(pElem=pRoot->FirstChildElement("WakeRelaxation")){
		

		TiXmlElement *pParamElement;
		
		if(pParamElement=pElem->FirstChildElement("ReplaceWake")){
			
			if (pEl=pParamElement->FirstChildElement("InitialWakeDirection")){
				WakeDir=ReadVector(pEl);
			}
			assert(pEl=pParamElement->FirstChildElement("StartSegmentLength"));
			StartSegmentLength=atof(pEl->GetText());
			
			assert(pEl=pParamElement->FirstChildElement("EndDistance"));
			EndDistance=atof(pEl->GetText());
			
			assert(pEl=pParamElement->FirstChildElement("NextElementRatio"));
			NextElementRatio=atof(pEl->GetText());
			
		}
		
		assert(pEl=pElem->FirstChildElement("NumberOfIterations"));
		NIter=atof(pEl->GetText());

		assert(pEl=pElem->FirstChildElement("RelaxationFactor"));
		RelaxationFactor=atof(pEl->GetText());

		assert(pEl=pElem->FirstChildElement("VortexCoreSize"));
		VortexCoreSize=atof(pEl->GetText());
		
	}
	
	
	
	assert(pElem=pRoot->FirstChildElement("Velocities"));

	for(pEl=pElem->FirstChildElement("Velocity"); pEl!=NULL; pEl=pEl->NextSiblingElement()) {

		
		
		Mesh=OrigMesh;

		Vector3D Velocity=Normalize(e1());
		
		Velocity=ReadVector(pEl);
		
		Mesh.ReplaceAllWakePolys(WakeDir, StartSegmentLength, EndDistance, NextElementRatio, WeldingTolerance);

		VortSystem System(Mesh);

		TiXmlElement *pPlanes;
		
		if (pPlanes=pRoot->FirstChildElement("SymmetryPlanes")) {
			TiXmlElement *pS;
			for(pS=pPlanes->FirstChildElement("Plane"); pS!=NULL; pS=pS->NextSiblingElement()) {
				Vector3D Normal,Position;
				Position=V0();
				Normal=e3();
				TiXmlElement *pParamElement;
				if(pParamElement=pS->FirstChildElement("Position")){
					Position=ReadVector(pParamElement);
				}
				if(pParamElement=pS->FirstChildElement("Normal")){
					Normal=ReadVector(pParamElement);
				}
				VortSymmetryPlane Plane(Position,Normal);
				System.AddSymmetry(Plane);
			}
		}

		System.CheckMesh(WeldingTolerance);
		
		System.SetVelocity(Velocity);
		
		
		
		for (int i=0; i<NIter; i++) {
			System.SolveSystem();
			
			if(i!=NIter-1){
				System.RelaxWake(RelaxationFactor,VortexCoreSize);		
			}
		}
		
		//	System.GetMesh().ExportMathematicaList("MathWing.dat");
		//	System.GetMesh().ExportCenterPoints("MathCenter.dat");
		
		/*	std::vector<Vector3D> Stream;
		 Vector3D R0=0.3*e2()-0.8*e1()-0.05*e3();
		 System.StreamLine(R0, 0.01, 10, Stream);	
		 WritePoints(Stream, "TestPoints.dat");
		 System.StreamLineRegular(R0, 0.01, 10, Stream, 0.01);	
		 WritePoints(Stream, "TestPointsRegular.dat");*/
		
		Memory Pressures,Gradients;
		MemoryVector Velocities;
		
		Vector3D F,M;
		
		System.CalculateDerivativePressures(Pressures,Velocities,&F,&M);
		
		//cout<<"Force= "<<F.v[0]<<", "<<F.v[1]<<", "<<F.v[2]<<",  Moment= "<<M.v[0]<<", "<<M.v[1]<<", "<<M.v[2]<<endl;
		cout<<"{"<<Velocity.v[2]/Velocity.v[0]<<", "<<F.v[0]<<", "<<F.v[2]<<"},"<<endl;
		
		string PLYName="finalmesh.ply";

		TiXmlElement *pSub;

		if (pSub=pEl->FirstChildElement("FinalMeshName")) {
			PLYName=pSub->GetText();
		}
		
		System.GetMesh().SavePLY(PLYName);
		
		System.CalculateDirectedAbsolutePressureGradient(Pressures,Velocities,Gradients);
		
		
		if (pSub=pEl->FirstChildElement("DipolesName")) {
			string PName=pSub->GetText();
			System.ExportData(PName,System.GetDipolesReference());
		}
		
		if (pSub=pEl->FirstChildElement("PressuresName")) {
			string PName=pSub->GetText();
			System.ExportData(PName,Pressures);
		}
		
		if (pSub=pEl->FirstChildElement("PressureGradientsName")) {
			string PName=pSub->GetText();
			System.ExportData(PName,Gradients);
		}
		
		//System.ExportData("TestPressures.dat",Pressures);
		//System.ExportData("TestGradients.dat",Gradients);
		//System.ExportIsWake("TestWake.dat");
		//System.ExportClC("TestClC.dat");
		
		//	cout<<"Saving="<<Mesh.ExportMathematicaList("TestWing.dat")<<endl;
		//	cout<<"BadEdgesExport="<<Mesh.ExportBadEdges("BadEdges.dat")<<endl;
	}

	
}
